Subspace representations in ab initio methods for strongly correlated systems 
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We present a generalized definition of subspace occupancy matrices in ab initio methods for 
strongly correlated materials, such as DFT+?7 and DFT+DMFT, which is appropriate to the case 
of nonorthogonal projector functions. By enforcing the tensorial consistency of all matrix opera- 
tions, we are led to a subspace projection operator for which the occupancy matrix is tensorial and 
accumulates only contributions which are local to the correlated subspace at hand. For DFT+J7 in 
particular, the resulting contributions to the potential and ionic forces are automatically Hermitian, 
without resort to symmetrization, and localized to their corresponding correlated subspace. The 
tensorial invariance of the occupancies, energies and ionic forces is preserved. We illustrate the effect 
of this formalism in a DFT+?7 study using self-consistently determined projectors. 
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I. INTRODUCTION 

The routine ab initio study of strongly correlated 
systems, that is those for which the accurate descrip- 
tion of the physics is beyond the capacity of mean-field 
methods such as Kohn-Sham density functional theory 1 
(DFT) within local or semi-local approximations to the 
exchange-correlation (XC) functional, remains a chal- 
lenge for electronic structure calculations. 

A number of sophisticated methods to correct the de- 
scription of strong correlation effects within DFT have 
been developed which provide a good compromise be- 
tween accuracy and computational expense. Successful 
examples include calculations using self-interaction cor- 
rected XC functionals 2 , exact exchange in DFT 3 and the 
GW approximation 4 , among others. Here we focus on 
methods, notably DFT + Hubbard U (DFT+L7) 5 and 
DFT + dynamical mean field theory (DFT+DMFT) 6 
for static and dynamical spatially-localized correlation 
effects, respectively, which share a common history and 
conceptual motivation based on models for Coulomb in- 
teractions such as the renowned Hubbard model 7 . 

In such methods, the electronic system is subdivided 
into a set of spatially-localized correlated subspaces, the 
description of the interactions in which is deemed to be 
beyond the capacity of the approximate XC functional, 
and the remainder which acts as a bath for particle ex- 
change and for which, due to its having a large kinetic en- 
ergy relative to Coulomb repulsion, the approximate XC 
functional performs adequately. In this manner, a model 
interaction may be used to augment and improve the de- 
scription of the screened Coulomb interactions between 
densities in the correlated subspaces while retaining the 
computationally inexpensive XC approximation for the 
remainder of the system. Generally for these methods, 
the occupancy matrix of each correlated subspace is the 
object which provides the necessary information on the 
electronic density to the model describing intra-subspace 
interactions. Denning the occupancy matrix of a corre- 



lated subspace using a set of orthonormal projectors is 
quite straightforward, yet the question of how to prop- 
erly extend the formalism to allow for the possibility of 
nonorthogonal spanning functions is one under active de- 
bate 8,9 and one of immediate practical consequence. 

It is frequently useful to permit the nonorthogonal- 
ity of the basis functions for the Kohn-Sham 1 states 
in ab initio methods which make use of sophisticated 
spatially-localized orbitals for such functions, particu- 
larly in linear-scaling density functional theory meth- 
ods 8,10 ~ 12 . Additionally, either for reasons of computa- 
tional convenience, as in Refs. 8, 13, and 14, or for the 
purposes of achieving self-consistency over the correlated 
subspaces, as in Ref. 15, it is common to use a subset 
of these nonorthogonal basis functions as projectors for 
the correlated subspaces, the subset termed Hubbard pro- 
jectors. We demonstrate here, however, that it may be 
hazardous to over-identify the Hubbard projectors with 
the basis set from which they are drawn. 

In this Article, we offer a revised definition of the sub- 
space occupancy matrix for ab initio methods which use 
nonorthogonal projectors to define the strongly corre- 
lated subspaces. We show that, by enforcing the ten- 
sorial consistency of all matrix operations, we are led 
immediately to a simple definition of the projection op- 
erator for each subspace which is fully localized to that 
subspace. In contrast to previously proposed formalisms 
of Ref. 8 and references therein, this gives rise to Hermi- 
tian corrections to the potential and ionic forces, without 
any post hoc symmetrization, which are also localized to 
the spaces in which the correlation correction is required. 
The resulting occupancy matrix reproduces the electron 
number of the subspaces and is tensorial. Thus, for ex- 
ample, its trace is invariant under both unitary rotations 
and the generalized Lowdin transformations 16 of Ref. 9. 

In order to illustrate the performance of the proposed 
formalism, we applied it to the DFT+t/ method in a 
study of two strongly correlated systems, namely bulk 
nickel oxide and the gas-phase copper phthalocyanine 
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dimer, with comparison to the most comprehensive al- 
ternative formalism available at the time of writing, the 
"dual representation" of Rcf. 8. A set of nonorthogo- 
nal generalized Wannier functions 10 , optimized using the 
projector self-consistent DFT+J7 method described in 
Rcf. 15, was used in order to carry out our computa- 
tional study with a minimum of user intervention in the 
construction of the nonorthogonal Hubbard projectors. 



II. NONORTHOGONAL REPRESENTATIONS 
OF THE OCCUPANCY MATRIX 

Generally, in order to extract low-energy Hubbard- 
model like models from ab initio DFT simulations, we re- 
quire the projection of the single-particle density-matrix 

/> w = £l^V£Wl' (i) 

where ipiu ^ s a Kohn-Sham eigenstate for spin channel 
a with band index i, crystal momentum k and occu- 
pancy f£\ onto a set of spatially localized subspaces. 
These subspaces where I is the site index, encom- 
pass that part of the Hilbert space of the Kohn-Sham 
orbitals which is deemed to be responsible for strong lo- 
calized Coulomb interactions beyond the scope of the ap- 
proximate XC functional. 

The occupancy of subspace which is delineated by 
a set of potentially nonorthogonal spanning projec- 
tors \^m)-, m G {1, . . . , M( J )}, dubbed Hubbard projec- 
tors, which are associated with subspace I, is generally 
given by the subspace-projected density matrix 

a(/)(<0 = p(/)t£(<0p(O ( 2 ) 

The Hubbard projection operator p(*\ the resolution of 
the identity for the space is defined in terms of the 
Hubbard projectors, but the exact manner in which this 
definition should be made has been the subject of some 
discussion, as we describe in the following. 

Some important conditions should be satisfied by a 
sound definition of the occupancy matrix of each corre- 
lated site, namely: all operations such as matrix products 
and traces should be tensorially consistent so that the 
total energy, potential and forces are tensorial invariants 
(unaltered by arbitrary transformations of the basis on 
which the projectors for that site arc defined); any poten- 
tial depending on that occupancy matrix should be Her- 
mitian and its action should be strictly localized to the 
correlated subspace while depending only on occupancies 
which are themselves localized to that subspace; the trace 
of the occupancy matrix should exactly reproduce the oc- 
cupancy of the correlated manifold on that site and if the 
site is extended to encompass the entire system then the 
total electron number should be obtained. 



A. The "full" and "on-site" representations 

We generally assume that a set of complex, mutually 
nonorthogonal Hubbard projectors are used for each in- 
dividual site and that the correlated subspaces possibly 
overlap (we do not consider transformations among the 
projectors of different correlated sites). Dual vectors of 
the Hubbard projectors must be defined with respect 
to some Hilbert superspace of the correlated manifold, 
%(!) z) C( J ), some possibilities for which are the subspace 
itself (i.e., T-L^ = C^), the union of all correlated sub- 
spaces (i.e., = {JjC^) and the space S spanned by 
all basis functions in the simulation cell (i.e., %W = S). 
The Hubbard projector duals are then generally given by 

\* (I)m )= E i*4 J) >s (7)aro > ( 3 ) 

aenc') 

where S^" is the contravariant metric tensor for the set 
of functions spanning %^ (the inverse of their overlap 
matrix). Physically meaningful inner products, e.g., ten- 
sorial invariants such as occupancies, energies or forces, 
are computed between functions and elements of their set 
of dual functions only (in the orthonormal case there is no 
practical distinction between functions and their duals). 
For a more detailed exposition of tensor calculus applied 
to problems in electronic structure theory, we refer the 
reader to Rcfs. 17. 

It is immediately clear that the simplest definition of 
the occupancy matrix for a given site, that is the projec- 
tion = J2mecc> \<Prn){ l Prn\ 01 the valence manifold 
over the site's Hubbard projectors, 

^ = (^1^)1^), (4) 

is invalid for nonorthogonal projectors. This widely-used 
definition of the occupancy matrix, which is entirely ap- 
propriate in the orthonormal case, such as calculations 
described in Ref. 18 and numerous citations therein, sim- 
ply neglects all nonorthogonality; the trace or powers of 
such a fully covariant tensor are not physically meaning- 
ful in the nonorthogonal case. 

A total site occupancy defined as a trace operation on 
this matrix, as in 

A^^E^' (5) 

m 

implies that such an occupancy is not, in general, a ten- 
sorial invariant since it is formed by a tensorially invalid 
summation over two covariant indices - as opposed to 
a meaningful contraction of indices of opposite tensor 
character. Occupancies, just like total energies, should 
be tensorial invariants, scalars which are unchanged by 
transformations of the basis on which the projector func- 
tions are defined. 

Progress was made in the definition of the occupan- 
cies of correlated subspaces via nonorthogonal projectors 
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when it was noted 14 that tensorially contravariant pro- 
jector duals should be involved, a concept known in other 
contexts for some time 19 . A definition of the occupancy 
matrix fully in terms of Hubbard projector duals was 
described in Ref. 14, for example, where the projection 
operator defined as = Y!,meC( I '>\'P^ I ' >m )('P^ m \> P r0 " 
vides an occupancy matrix 

= S^ ma (ip^\p^\ipf)S^ 0m '. (6) 

The indices a and j3 run over the spanning vectors of 
the contravariant metric (i.e., the inverse overlap matrix) 
on a superspace T~L^ of the correlated manifold 
C^. Here and hereafter, we make use of the summation 
convention 20 , whereby repeated indices within the same 
expression are summed over unless in parentheses. 

Unfortunately, the matrix trace and powers of Eq. 6 are 
not tensorially valid, as can be seen by taking the exam- 
ple of the square of this contravariant occupancy matrix, 
which is of interest for density-density self-interaction 
corrections to approximate XC functionals. The result- 
ing expression for the squared occupancy matrix 



n 2{I)(a)mm' = ^ n (I)(a)mm" n (I)(a)m" m' ^ 
m"eC( J ) 

implies that the operator 



p (I) = E w {I)m ")(v {I)r ~" 

m"£CV) 



(8) 



forms a tensorially traceable identity on It does not 
in the case of nonorthogonal projectors, however, since 
an identity operator may only be formed via the outer 
product between a projector and a projector dual, and 
not a dual vector and its own complex conjugate, includ- 
ing the case where the correlated subspace is extended to 
the Hilbert space of the entire system (C^ = S). 

The shortcomings in the two definitions of the occu- 
pancy matrix described above have been previously de- 
scribed in detail by Han ct. al. in Ref. 8 and are dubbed, 
respectively, the "full" (Eq. 4) and "on-site" (Eq. 6) 
representations in the nomenclature described therein. 
These authors concentrated on the special case where the 
dual-generating superspace H {1) is the space spanned by 
all basis functions {|</> Q }} in the simulation cell, so that 
U {1) = S, in which case S a f) — (4> a \4>i3) an d the Hubbard 
projectors form a subset of the basis set. Thus, the same 
contravariant metric for all basis functions in the simu- 
lation cell is used to generate the dual functions on each 
correlated site and in this case the "full" and "on-site" 
occupancy matrices simplify, respectively, to 



and 



WW = E S mecWa K^S [im , eC{1) (9) 
a, pes 



Here, K^ al3 — [<f) a \p^\(j)^) is the representation of the 
density matrix in terms of basis-set duals, known as the 
density kernel. The notation S meC (i) a reminds us that 
m and a run over the spanning vectors of two different 
spaces, and = S, respectively, so that the block 
of S„ in question is generally not square. 



B. The "dual" representation 

Han et. al. 8 , whose invaluable contribution on this 
subject addressed many of the salient issues, pointed out 
that the total number of electrons is not recovered by 
the trace of the occupancy matrix if the site is extended 
to include the entire simulation cell using the "full" and 
"on-site" representations. They proposed an alternative 
"dual" representation which solves this particular prob- 
lem and is generated by the projector 



p (I) = - 2 E (i^x^i + i^)^ 



and the corresponding occupancy matrix 



(11) 



A E 



jy-(cr)m£C<- I} a q 



(12) 



In the "dual" representation, the contravariant metric 
on the basis set is used to form the Hubbard projector 
duals (which are therefore delocalized across the entire 
simulation cell, in general, since the inverse overlap ma- 
trix is dense even when the overlap matrix itself is sparse) 
via \ip^ m ) = ^ ae g\tpa^) S am . Symmetrization is then 
carried out in order to both provide a symmetric occu- 
pancy matrix and to recover a Hermitian potential. 

The "dual" representation shares with the "full" repre- 
sentation the attribute of Hcrmiticity and, furthermore, 
it has a tensorially and physically meaningful trace. As 
such, to our knowledge, it provides the most favourable 
occupancy definition hitherto available. However, this 
occupancy matrix is tensorially ambiguous, consisting of 
the sum of tensors of differing index character. 

One cannot generally symmetrize or antisymmetrize a 
tensor over indices of mixed covariant and contravariant 
character in this way and obtain a matrix which trans- 
forms as a tensor (one which may be used to generate 
tensorially invariant occupancies, local moments or ener- 
gies). Thus, while providing a significant improvement 
over previously suggested definitions of the occupancy 
matrix due to its tensorially invariant trace, the "dual" 
representation suffers similar problems with matrix pow- 
ers as other representations: if we attempt to compute 
the square of this matrix we obtain tensorially inconsis- 
tent, and thus physically meaningless, terms in the prod- 
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C. Requirement for a subspace-localized Hermitian 
projection operator 

Let us step back for a moment and consider why any 
projection operator of the form 

P (I) = E \^ ] )S am (^\ (13) 
mec<'> 

requires symmetrization to the "dual" form in order to 
provide a Hermitian potential operator. 

An arbitrary potential operator V, operating on 
the subspace which could represent the screened 

Coulomb interaction, for example, has matrix elements 
in the frame of Hubbard projectors given by 

V^ m ' = J2 (V$\V\^)S™' . (14) 

The potential operator is easily shown to be non- 
Hermitian in the case where m,m' £ C %^ C S, 
and ^ strictly holds, since 

V(i) = pwvp(i) = lip g) )v (i) mm ' {ip (i) tl (15) 

= E \^)s (I)ma v$sw™\^,\ 

= w (I)m )v^l^ {I)m '\ = pmpw = v^. 

The reason for this non-Hcrmiticity is that the indices 
a, f3 do not generally run over functions spanning just the 
correlated space but rather over those that span a 
superspace e.g., typically over the basis functions 

in the simulation cell, %W = S. This observation is 
quite general: the dual projectors must be constructed 
using the metric on precisely the space spanned by the 
projectors themselves in order to build a Hermitian pro- 
jection operator and hence a Hermitian potential. This 
cannot be circumvented in a tensorially-consistent way 
by symmetrizing operators since tensors can only be sym- 
metrized over pairs of indices if they are either both of 
covariant character or both of contravariant character. 



D. The "tensorial" representation 

Based on the above, we deduce that, in order to build a 
tensorially-consistent occupancy matrix which generates 
a Hermitian potential, the projection operator for a given 
subspace must necessarily be constructed using exact 
dual Hubbard projectors with respect to that subspace 
only. Thus, with the covariant overlap matrix of Hubbard 
projectors defined by 

0% m , = {^\^,), (16) 



that is an individual x covariant metric tensor 
for each correlated site /, the proper dual vectors \if^ m ) 
are constructed using the corresponding contravariant 
metric 0( 7 ) m m , as per 

\^ m )= E \^)0 (I)m ' m . (17) 

Here, we emphasize that the contravariant metric is ob- 
tained via a separate x inverse operation for 
each site, so that 

OW m ' m "0$, m = Si n . (18) 

In the special case where the Hubbard projectors are 
drawn from the set of functions used to represent the 
Kohn-Sham wave-functions, the overlap matrix of duals 
f or each site cannot generally be extracted im- 
mediately from the metric S" on S. However, in this 
particular case, the 0[ ! J matrix for each site is merely 
a sub-block of the basis-function overlap S», and, from 
this, the contravariant O^** for each site can be com- 
puted by a separate inverse operation for each site which 
is typically fast, due to the small matrix dimension. 

Employing this definition of the metric tensor on each 
subspace, the projector duals remain manifestly as local- 
ized to the correlated subspace as the projectors them- 
selves, they pick up only subspace-localized contributions 
to the occupancy and can only apply subspace-localized 
corrective potentials, as we would expect for local self- 
interaction corrections such as DFT+J7 or its extensions. 
The Hubbard projection operator, in what we will denote 
the "tensorial" representation, 

P (I) = E \<P%)O m ' m {<P$\ (19) 

is Hermitian and thus gives rise to a Hermitian potential, 
without resort to symmetrization, since is a square 

overlap matrix, 

V(D = pWvPW = \^)V^ mm '(^}\ (20) 

= \^W I)mm 'v^l^ I)m '' m ''\^l\ 
= W {I)m )v^Av {I)m '\ = P^vpw = vw. 

The occupancy matrix is most easily expressed in its 
singly covariant and singly contravariant form, though 
other forms are readily obtainable from the metric tensor, 
so manipulations of the following form can be made: 

n,' = 0..n" = n..O" = 0..n',0". (21) 

The contravariant-covariant or covariant-contravariant 
forms of the tensorial occupancy matrix, the latter given 
by (the second line applies to the special case where Hub- 
bard projectors are drawn from the basis set) 

= E S ma K^S (im n0^ m " m \ (22) 
a, pes 
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possess a common tensorially invariant trace (a proper 
contraction over one covariant and one contravariant in- 
dex) which recovers the exact number of electrons in the 
correlated subspace by construction (the so-called sum- 
rule), so that 



E 



E '< 

neC'> 



(23) 



Their powers themselves remain tensors, for example 

the square n 2 ^^™" 1 = n^ I ^ rT m m n^^,T is itself a well 
behaved singly covariant and singly contravariant tensor 
with an invariant trace. This occupancy matrix trace is 
easily demonstrated to be invariant under rotations of the 
set of Hubbard projectors on its site and it is independent 
of the basis used to represent the Kohn-Sham states. 

An occupancy matrix which is invariant under element- 
wise transpose might lend itself to an interpretation as 
quantifying the charge shared between Hubbard projec- 
tors, and indeed it does in the case of orthonormal Hub- 
bard projectors. However, it is worth emphasizing that 
in the case of a set of nonorthogonal Hubbard projectors, 
these functions are merely spanning vectors with no rig- 
orous physical meaning and, generally, no such interpre- 
tation of charge shared between orbitals may be safely 
made. In fact, in the nonorthogonal case, the occupancy 
matrix should not be generally expected to be invariant 
under element- wise transpose, i.e., n^ 1 ^ n^ m = n m m . 
Rather, if the duals are defined in a way which pre- 
serves the tensorial consistency of inner products, the 
occupancy matrix must satisfy instead the more general 
expression = O mm »n m m ,„O m m , where O.. and 

O" are the covariant and contravariant metric tensors, 
respectively, on the subspace in question. As a result, 
only the diagonal elements of the occupancy matrix can 
be imbued with a intuitive meaning in the sense of oc- 
cupancy; symmetrizing the matrix does not recover such 
an interpretation for the off-diagonal elements. 



transition-metal or Lanthanoid atom, is used to define 
the occupancy matrix of the correlated subspace at each 
site. The particles occupying these subspaces interact 
strongly with each other by comparison with their inter- 
action with the bath; each subspace acts as an open quan- 
tum system. As such, we may separately impose the Fock 
antisymmetry condition for the projected wave-function 
for each strongly correlated subspace, so that the sub- 
space occupancy, the projection of the full single-particle 
density onto the subspace in question, should itself be a 
valid density-matrix operator (it should be idempotent 
and reproduce the electron number of that subspace). 

Since the idempotency of the density matrix for the full 
system is a condition which must be exactly satisfied, and 
the idempotency of correlated sites is a competing con- 
dition (the Hubbard projectors differ from Kohn-Sham 
orbitals in general), the subspace idempotency may be 
only partially enforced, for each correlated site, using an 
idempotency penalty functional of the form 



la 



n (/)(<x)_a(j)(<02> 



(24) 



which disfavours deviation from wave-function anti- 
symmetry in the strongly correlated subspaces. The 
pre-multiplier \W< a ^ is usually approximated by a sin- 
gle scalar for each site, where it is identified as 



(25) 



half of the screened, subspace-averaged Coulomb interac- 
tion. If we further assume an orthonormal set of Hubbard 
projectors for each site, the functional is easily recognis- 
able as the familiar rotationally invariant DFT+U cor- 
rection term of Cococcioni and de Gironcoli in Ref. 18, 



E 

la 



n 



mm' "' in' rn 



(26) 



III. APPLICATION TO THE DFT+C/ METHOD 

In this section, we illustrate the practical application 
of the "tensorial" representation to a particular method 
for strongly correlated materials, namely the simplified 
rotationally invariant DFT+U correction of Refs. 18 and 
21. We provide the necessary expressions for the tensori- 
ally invariant DFT+U terms in the energy, potential and 
ionic forces for use with nonorthogonal Hubbard projec- 
tors, which is of some importance since such a set is often 
used in contemporary high-accuracy, particularly linear- 
scaling, implementations 8 ' 10-12 . We have recently shown 
that an efficient set of Hubbard projectors can be con- 
structed which is self-consistent with the set of truncated 
nonorthogonal generalized Wannicr functions which min- 
imize the DFT+U total energy 15 . 

In the DFT+U approach, a set Hubbard pro- 

jectors, typically spatially localized on a particular 



A. The tensorially invariant DFT+(7 functional 

Let us consider how we might generalise this DFT+U 
penalty-functional to accommodate an orbital-dependent 
interaction tensor. The Coulomb interaction tensor U for 
a given spin channel and site (considering the same Hub- 
bard projectors for different spins for brevity of notation) 
is given generally by the two-centre integral (N.B. using 
the Dirac, and not Mullikcn, convention) 

U ( ±„ = <«^ (7)(CT) M \^WS). (27) 

Here, U^ (r, r') is the Coulomb interaction screened ac- 
cording to mechanisms described by an appropriate the- 
ory such as linear-response 13,18 ' 21 , constrained LDA 22 , 
constrained RPA 23 or constrained adiabatic LDA 24 . 
Coulomb repulsion is represented by those terms for 
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which m = m"; m! = m'", while direct exchange is given 
by those elements with m = m'";m? = m". 

In the general, nonorthogonal case, care must be taken 
in employing the U tensor in order to preserve the ten- 
sorial invariance of the DFT+[7 energy. For example, 
if a tensorial invariant is required which provides the 
sum of the part of the tensor describing density-density 
Coulomb repulsions, it should correctly be computed by 
contracting covariant and contravariant indices in pairs 
of indices of opposite character, i.e., double-sums of the 
form, where m, ml <G {!,..., M^}, 



U, 



mm 
mm' 



u 



mm 
mm' i 



u, 



m m 
% m' ) 



u 1 



(28) 



are admissible while those of the form U mm i mm i or 

jjmm mm break tcnsoria l 

invariance. Indices are raised 
and lowered simply using the metric tensor of the cor- 
related subspace to which the U tensor corresponds, the 
contravariant O" or covariant O.., respectively, e.g., 



U„ 



= O" 



U, 



mm" m' m 



(29) 



Purely as an illustration of this principle, a simple 
projector-decomposed tensorially invariant penalty func- 
tional may be constructed using pairwise contractions 
over the four indices, as in 



E 

la 



1 mm' 

2 mm" 



*X m _ 

"ml" Ur. 



n„ 



(30) 



A commonly used approximation for the screened 
Coulomb interaction, that which we use, is where the 
interaction tensor (itself an inverse response function) is 
averaged over the subspace (i.e., both over perturbing 
and probing indices), providing a scalar density-density 
Coulomb interaction. The usual DFT+l/ penalty func- 
tional in this fully averaged approximation is thus given, 
in tensorially-invariant form, by the expression 



E 

la 



1 uW 7 

2M( 7 ) 2 m m 



(31) 



B. DFT+[7 potential and ionic forces in the 
tensorial formalism 



The DFT+[7 term in the Kohn-Sham potential, gen- 
erally given (for real valued U tensors) by 

=J2\^ I)m )Vi I)Mm '(v { J\, (32) 



has matrix elements, in the case of averaged U, given by 



y(I)(<r)m' _ \ jA 1 ) ™"m" 
m 2_M(-f) 2 m " m "' 



5™' - 2n^ m ' 



The DFT+U potential is Hermitian by construction 
when the Hubbard projection operator built with the 
subspace-local tensorial representation of Eq. 19, is used. 



No symmetrization of the occupancy matrices is then 
needed to ensure this Hermiticity and the potential acts 
strictly within the spatial extent of the subspace of whose 
occupancy it depends. 

Correspondingly, the DFT+U contribution to the force 
on the ion labelled J, with position Rj, is given by 



^ = -E 

la 

-E 



/ dtp 



(I) 



\ dRj 



JI)(a)m'/ (I) 



ln (I)\ n (I)m'm" (I)(a)m'" ,AI)(a)r 
fm' / U U m" V m"> 

(I) 



m" \ri{I)m"m"' *AI)(cr)m 

, )v v m »> ■ 



Here, we have made simplifications such as 



dip 



dR 



m' \/~\{I)m' m 

It I 



(33) 



which is valid if the subspace metric tensor is position 
independent, in particular if the Hubbard projectors are 
simply spatially translated when their host ion is moved. 
Our force expression holds, of course, only if we are 
on the Hellmann-Feynman surface, where the density- 
matrix commutes with the Hamiltonian. 



IV. BULK NICKEL OXIDE 

The first row transition-metal monoxide NiO poses 
some difficulties to Kohn-Sham density-functional theory 
and to electronic structure theories generally. As such, it 
has served as a valuable proving-ground for approaches 
such as periodic unrestricted Hartree-Fock theory 25 , the 
self-interaction corrected local density approximation 2 , 
the GW approximation 4 and DFT+DMFT 26 . Exper- 
imentally, the paramagnetic phase of NiO is found to 
possess a rock-salt crystal structure with a lattice con- 
stant of approximately 4.17 A 27 . At ambient tempera- 
ture, NiO is a type-II antiferromagnetic insulator with a 
local magnetic moment of between 1.64 fis and 1.9 ^s 18 - 
Due to the persistence of its magnetic moment and op- 
tical gap, which is approximately 4 eV, above the Neel 
temperature, it falls broadly into the category of a Mott 
insulator 25 with a charge-transfer insulating gap of pre- 
dominantly oxygen 2p to nickel 3c?-orbital character 25 ' 28 . 
It has long been recognised that LDA-type exchange cor- 
relation functionals 29 qualitatively fail to reproduce the 
physics of this material, grossly under-estimating the lo- 
cal magnetic moment, the Kohn-Sham gap and assign- 
ing an incorrect fully nickel 3(i-orbital character to the 
valence band edge. We stress, however, that the Kohn- 
Sham gap is not comparable to the experimental insulat- 
ing excitation gap, even for the exact XC functional 30 . 

The DFT+t/ method has previously been applied, in 
numerous incarnations, to bulk NiO and it is known to 
recover the principal features of this strongly-correlated 
oxide 5,13 ' 18,27 ' 31 . Moreover, generalizations to DFT+U 



such as first-principles methods for calculating the Hub- 
bard U parameter 13 ' 18 , the BFT+U+V method for in- 
cluding inter-site interactions 32 and, most pertinent for 
this study, previous investigations into subspace rep- 
resentations of nonorthogonal Hubbard projectors in 
DFT+C/ 8,9 have also been applied successfully to this sys- 
tem. We have chosen to study NiO, therefore, because 
it is so well characterized and we have performed cal- 
culations which we hope will be complementary to those 
described in Ref. 8, where the "full" , "on-site" and "dual" 
representations of a linear-combination of pseudo-atomic 
orbital basis were compared. 

A. Computational methodology 

Calculations of the ground-state electronic structure 
of bulk antifcrromagnetic nickel oxide were carried out 
within collinear spin-polarized Kohn-Sham DFT 1 , and 
the simplified DFT+U method 18 . The linear-scaling 
ONETEP first-principles package, described in detail 
in Refs. 33, was used. The LSDA (PZ81) exchange- 
correlation functional 29 , with norm-conserving pseu- 
dopotentials 34,35 , was invoked throughout. Periodic 
boundary conditions were used with a 512-atom super- 
cell and the Brillouin zone was sampled at the T-point 
only. A systematic variational basis of Fourier-Lagrange, 
also known as periodic cardinal sine or psinc, functions 36 , 
was used, equivalent to a set of plane- waves bandwidth 
limited to a kinetic-energy cutoff of 825 eV. 

In the ONETEP method, the Kohn-Sham density- 
matrix is represented in the separable form 

p^(r,r') = (Mr)X (CT)a %(r') (34) 

in terms of a set of covariant nonorthogonal general- 
ized Wannier functions (NGWFs) 10 , {0. (r)}, and a cor- 
responding contravariant density kernel, K" , for each 
spin channel. The density kernel was not truncated 
in the calculations described here. In the ONETEP 
method 33 , the total energy is iteratively minimized both 
with respect to the elements of the density kernel for a 
given set of NGWFs 37 , using a combination of penalty- 
functional 38 and LNV 39 techniques to ensure the validity 
of the density matrix, and with respect to the expansion 
coefficients of the NGWFs in the psinc basis. The con- 
verged NGWFs (a minimal set of nine functions for nickel 
4s, 4p and 3d and four for oxygen 2s and 2p, truncated 
to an atom-centered sphere of 4.0 A, were employed in 
calculations on NiO) are those which are optimized to 
minimize the total energy and are thus adapted to the 
chemical environment, incorporating all valence-electron 
hybridization effects in the ground-state density. 

Our principal purpose was to provide an appraisal of 
the difference in predicted electronic properties, if any, 
given by DFT+U when using nonorthogonal Hubbard 
projectors with either the "dual" or "tensorial" represen- 
tations of the correlated subspaces. The "dual" represen- 
tation, in particular, was selected for comparison since it 
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FIG. 1. (Color online) Isosurfaces of the set of nonorthogonal 
generalized Wannier functions (NGWFs) on a nickel atom 
in NiO. The NGWFs are those computed at projector self- 
consistency in the "tensorial" representation at LDA+U = 
6 eV. Those in the left column (predominantly 3d — tig char- 
acter) and the top and bottom NGWFs in the middle column 
(predominantly 3d — e g character) are those used as Hubbard 
projectors, while the remaining NGWFs (pseudized 4s-like in 
the centre and pseudized 4p-like in the right column) lie out- 
side the correlated subspace on that atom. The isosurface is 
set to half of the maximum for the 4s and 4p-like NGWFs 
and 10~ 3 times the maximum for the 3d-like NGWFs. 



appears to be the most sophisticated of the previously 
proposed subspace definitions - it has a tensorially in- 
variant occupancy matrix trace which cannot be said of 
the manifestly incomplete "on-site" and "full" represen- 
tations. The latter three representations were previously 
compared in detail in Ref. 8. 

The underestimation of the NiO lattice parameter with 
respect to experiment by the LDA, as well as the ability 
of DFT+U to correct this, for a particular U value, has 
been known for some time 27 . The U parameter required 
to correct the lattice will depend on details of the under- 
lying XC functional, the precise DFT+C/ functional used, 
the pseudopotentials and both the form of the Hubbard 
projectors and the definition of the subspace projection 
operators. In order to provide an unbiased analysis of 
different subspace definitions, therefore, we employed the 
same experimental lattice constant for all calculations. In 
order to obviate intervention in the construction of the 
correlated subspaces, so far as possible, we carried out the 
DFT+U calculations in the projector self-consistent for- 
malism described in Ref. 15. We also include, for the pur- 
poses of comparison, the results of conventional DFT+U 
calculations using hydrogenic 3e?-orbital Hubbard projec- 
tors (in which case there is no ambiguity in the represen- 
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FIG. 2. (Color online) The total occupancy of a correlated 
subspace in NiO (left axis) and the maximum element on the 
diagonal of the occupancy matrix (right axis) , as a function of 
the interaction U. Values are computed with orthonormal hy- 
drogenic Hubbard projectors and self-consistent NGWF pro- 
jectors in both the "dual" and "tensorial" representations. 

tation for a given choice of projectors) which were used as 
the initial guess for the projector self-consistency cycle 40 . 

In the projector self-consistent DFT+U scheme of 
Ref. 15, the set of five converged NGWFs of maximal 
3(i-orbital character on a transition-metal atom respon- 
sible for strong correlation effects are selected as Hub- 
bard projectors to redefine the DFT+U occupancy ma- 
trices for the total energy minimization in the next pro- 
jector iteration. The energy is not directly minimized 
with respect to the expansion coefficients of the Hub- 
bard projectors (since it would violate the variational 
principle if either the Hubbard projectors or U were al- 
lowed to change during energy minimization ), but the 
projectors are updated in a manner reminiscent of the 
density-mixing method for solving non-linear systems 41 , 
converging towards those which equal a subset of the 
NGWFs which minimize the DFT+U energy functional 
which they themselves define. The projector-update pro- 
cess alternates between direct variational minimization 
of the total energy and projector renewal until both are 
individually converged. 



B. Occupancies and magnetic dipole moments 

In agreement with a number of previous studies 9 ' 27,31 , 
we find that the LDA correctly favours antiferromagnetic 
ordering in NiO, albeit with diminished local magnetic 
moments and a greatly underestimated Kohn-Sham gap. 
The DFT+U correction enhances the antiferromagnetic 
order with increasing U , monotonically increasing the 
magnetic dipole moments. Also in accordance with previ- 
ous work 8,9 , we have found that the DFT+U occupancy 
matrix and local magnetic dipole moment associated with 
the correlated subspaces depend significantly on the def- 



FIG. 3. (Color online) The projection of the magnetic dipole 
moment onto the DFT+U correlated subspace on nickel 
atoms in NiO, computed as a function of the interaction U. 
Values are computed as in Fig. 2. 

inition of the correlated subspace projection operator. 

Turning first to the total occupancy of the correlated 
subspaces, shown in Fig. 2, we find a steady decrease with 
increasing U parameter, which is almost entirely due to 
the DFT+U correction introducing a repulsive potential 
to the less-than-half occupied nickel 3d— e g orbitals of the 
minority spin channel. Conversely, we notice that for the 
largest element on the diagonal of the occupancy matrix 
(which is almost identical to that of the other orbitals of 
the same symmetry), DFT+U introduces an attractive 
potential that tends to fully occupy the corresponding 
orbital. 

The maximal occupancy element for hydrogenic pro- 
jectors, for those projectors most commonly used in 
DFT+U which arc not adapted to their chemical envi- 
ronment and so cannot fully account for densities devi- 
ating from spherical symmetry, slowly approaches unity 
and we conjecture that a rather excessive U value would 
be needed to complete the orbital filling. On the other 
hand, if we look at self-consistent NGWF projectors in 
the "dual" representation, there is a tendency to overfill 
the most fully occupied Hubbard projectors, to wit the 
occupancy begins to exceed unity beyond U ~ 3 eV. This 
latter affliction is a rather hazardous one for the DFT+U 
functional, since the contribution to the energy correction 
arising from orbitals exhibiting it may become negative 
in severe cases; this is incorrect behaviour for a penalty- 
functional in any case. The reason behind this excessive 
occupancy is the spurious non-locality of the Hubbard 
projector duals in the "dual" representation, they may 
pick up density contributions from all across the simula- 
tion cell. On the contrary, when self-consistent projectors 
are used in the "tensorial" representation, the maximal 
matrix elements tend asymptotically to unity with in- 
creasing U, as expected (reaching 0.9998 at U = 8 eV). 

In order to test the dependence of our computed DFT 
properties on the XC functional used for pseudopotential 
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generation, these dependencies known to be potentially 
substantial when non- linear core corrections are used 42 , 
we also performed our "hydrogenic" calculations using 
LDA 29 pseudopotentials with parameters closely match- 
ing to the GGA set 35 . Moving from the latter to the 
former pseudopotentials, we observed a reduction of the 
local magnetic moment by 0.007 {1b at U = eV, up 
to 0.02 he at U = 4 eV, whereat the reduction remains 
with further increase in U. The total correlated subspace 
occupancy is rigidly increased by approximately 0.02 e. 
The maximum occupancy matrix changes by no more 
than 0.001 electrons any U tested. The occupancy ma- 
trices thus depend on the choice of pseudopotential, as 
expected, as do derived properties, but not sufficiently to 
influence our observed trends, indeed by a small amount 
compared to the dependence on the U parameter and 
subspace projection definition. 

Considering the local magnetic moment on the nickel 
atoms, depicted in Fig. 3 and defined as the difference 
of the traces of the DFT+J7 occupancy matrices of the 
two spin-channels, we observe the expected increase with 
the U parameter as the majority and minority channels 
of the magnetization-carrying orbitals become increas- 
ingly filled or emptied, respectively. The NGWF pro- 
jectors in the "dual" representation yield greater local 
magnetic moments than the representation-independent 
hydrogenic projectors and, in turn, those are larger than 
the moments in the "tensorial" representation. Conse- 
quently, we would expect the exchange splitting which 
makes up a large contribution to the insulating gap in this 
material (it is well-described within unrestricted Hartree- 
Fock theory 25 ), to follow the same trend. While this 
behaviour may seem a somewhat unfavourable reflection 
on the "tensorial" representation, it is fully in line with 
our understanding that the "dual" representation (or 
any related delocalized "Mulliken" -type analysis) picks 
up additional contributions from magnetization densities 
of neighbouring atoms by construction. The previously 
demonstrated strong dependence of the moments on the 
definition of the subspace occupancy matrices 8 , taking 
the 4 eV spread of U values which approximately yield 
1.48 fiB in our calculations as an example, demonstrates 
the hazard incurred by comparing U parameters used 
with DFT+J7 methods that differ in their technical de- 
tails. 



C. Kohn-Sham eigenspectra 

The Kohn-Sham eigenspectrum computed for 
NiO using DFT+U with both our "best guess" 
system-independent hydrogenic projectors and self- 
consistently determined NGWF projectors agree closely. 
Moreover, in agreement with previous studies of the 
dependence on the occupancy matrix definition when 
using nonorthogonal Hubbard projectors 8 ' 9 , the repre- 
sentation dependence of spectral features is rather subtle 
and is considerably less significant than the dependence 




(Energy - E Fermi ) (eV) 



FIG. 4. (Color online) Top: Density of Kohn-Sham states per 
spin per atom of NiO at LDA+f/ = 6 eV, together with its 
projection onto the union of correlated subspaces using hydro- 
genic Hubbard projectors and NGWF projectors in the "ten- 
sorial" and "dual" representations. Bottom: The decomposi- 
tion, in the NGWF- "tensorial" representation, of the density 
of states for a given spin channel into its contributions from 
NGWFs on nickel atoms with magnetization aligned (major- 
ity) and anti-aligned (minority) spins, the correlated subspace 
projections of each, and the contribution due to NGWFs on 
oxygen atoms. 



on the U parameter. That is not to say, however, that 
the differences yielded may be guaranteed to be fully 
recovered by a self-consistent determination or arbitrary 
variation of the interaction U , since we observe different 
dependences on this parameter for different spectral 
peaks depending on the subspace representation. 

Considering, for example, a Hubbard U value within 
the range of values known to give reasonable agreement 
with experiment, namely U = 6 eV, we have shown the 
total density of states (DoS), and its correlated subspace 
projection, in the three representations of interest in the 
top panel of Fig. 4. The bottom panel shows the de- 
composition of the "tensorial" DoS into its contributions 
from oxygen atoms and both predominantly spin-aligned 
(majority) or spin-antialigned (minority) nickel atoms. 
Although all of the dominant features are shared between 
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FIG. 5. (Color online) The Hubbard U dependence of the 
Kohn-Sham band-gap of NiO at LDA+U. 



the eigenspectra of the various representations, there are 
some discrepancies which are worth noting. Most no- 
table is the trend for the insulating gap to open slightly, 
predominantly at the minority e g peak at « 2 eV, as 
we go from "tensorial" NGWF (2.35 eV) to hydrogenic 
(2.60 eV) to "dual" NGWF representations (2.68 eV). 
We attribute this to changes in the exchange splitting 
provided by the enhancement of the magnetic moment, 
which follows the same trend, as can be seen in Fig. 3. 
The localized character of the valence band edge is not 
significantly representation-dependent. 

We show the [/-dependence of the Kohn-Sham insu- 
lating gap in the three DFT+C/ correlated subspace def- 
initions in Fig. 5. In all cases, we recover the canoni- 
cal DFT+C/ description of this material. With increas- 
ing interaction parameter U the tendency is that: the 
low-energy (primarily majority-channel e g -like) peak falls 
deeper into the valence band as an attractive potential 
is applied to fill it completely; the strongly nickel ti g - 
like valence-band edge at the LDA level gives way to 
hybridized oxygen 2p character as the i2g-like states are 
pushed to lower energies; and the minority-channel Nickel 
e g -like first peak in the conduction band is increased in 
energy as its partial occupancy causes it to be subjected 
to a repulsive corrective potential. 

Overall, we reiterate that the effects on the spectra 
due to the local or non-local construction of the Hub- 
bard projector duals, at least in this material, are not 
sufficiently great to reasonably draw conclusions regard- 
ing the relative merit of methods based on agreement, 
or otherwise, with experimental observations. Rather, in 
this matter, points of principle such as the preservation 
of tensorial invariance, or the avoidance of occupancies 
exceeding unity (as observed in Fig. 2), must therefore 
take precedence in our view. 



V. COPPER PHTHALOCYANINE DIMER 

Open-shell molecular systems containing transition 
metal ions sometimes pose a challenge to first principles 
simulation within LDA-based approximations 43 . This 
is partially due to the tendency of such approximate 
XC functionals to excessively delocalize magnetization- 
carrying orbitals in such systems. As noted in Refs. 21, 
44-46, both energetic properties, such as magnetic cou- 
pling, and also spectroscopic features, such as the nature 
of the insulating gap and multiplet splittings, can con- 
sequently be poorly reproduced by such functionals. So- 
phisticated ab initio techniques such as the GW approx- 
imation and local correlation methods such as DFT+C/, 
whose traditional realm of application lies in extended 
systems such as extended oxides and their interfaces, are 
being increasingly applied to molecular systems and clus- 
ters (see for example Refs. 21, 44-47). 

It is thus of some importance, and perhaps timely, 
to consider molecular systems on a similar footing to 
solids when considering the merit of projection methods 
for DFT+C/. The correlated orbitals in molecular sys- 
tems may be rather more spatially diffuse and deviate 
further from spherical symmetry than their counterparts 
in solids. As a result, the issue of Hubbard projector- 
dependence in DFT+C/ and then the manner in which 
the projection operator is constructed from those projec- 
tors, particularly the degree of non-locality in the Hub- 
bard projector duals, can be expected to play a more 
significant role in the description of molecular systems. 

With a view to analyzing the dependence on the cor- 
related subspace definition, or occupancy representation, 
in the case of molecular systems, we applied our method- 
ology to the ground-state of a binuclear open-shell (an- 
tiferromagnetically coupled) singlet complex, the copper 
phthalocyanine dimer denoted a-Cu(II)Pc2. Crystalline 
CuPc is a semiconducting blue dye which, in pure thin- 
film form and more exotic derivatives, is currently at- 
tracting intense experimental and theoretical interest due 
to its potential for use as a flexible organometallic pho- 
tovoltaic material 48 , as part of field-effect transistors 49 
and, due to its magnetic functionality, in spintronic data 
storage or processing devices 50 . In this system, two cor- 
related subspaces delineated by copper 3c?-like states are 
spatially well separated, with approximately 3.77 A be- 
tween centres, and there is minimal electronic bonding 
between the localized orbitals in the open Cu-3d shells in 
the two approximately planar moeities. The result is a 
very weak indirect-exchange, i.e., acting via intermediary 
delocalized ligand states, S = h antiferromagnet with a 
Heisenberg exchange coupling constant of J w — 1.50K; 
for a detailed analysis of this mechanism see Ref. 51. 

It is important to emphasise that although correc- 
tive techniques for localized correlation effects such as 
DFT+C/ have been shown to be somewhat beneficial in 
the context of organometallic molecules 52,53 , the are by 
no means the only, or perhaps favourable, methods for 
such systems. For Cu(II)Pc2, as we go on to show, 
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FIG. 6. (Color online) Spin-density isosurfaces at 5% of max- 
imum in Cu(II)Pc2 at projector self-consistent GGA+U = 
6 eV within the NGWF- "tensorial" subspace representation. 



the magnetization-carrying copper 3c? orbitals are partly 
delocalised and thus not fully recovered by DFT+/7. 
Systems of this type have been described with partic- 
ular success, notably 51,54,55 , using hybrid XC functionals 
comprising a fraction of non-local Hartree-Fock exchange 
more appropriate to these molecules. 

The ONETEP method 33 was used, as before, with 
r-point Brillouin zone sampling, norm-conserving pseu- 
dopotentials 34,35 and a set of nine NGWFs (4s, 4p and 
3d) for copper ions, four each for carbon and nitrogen 
(2s and 2p) and one for hydrogen (Is). An NGWF 
cutoff radius of 5.3 A and an equivalent kinetic-energy 
cutoff of 1000 eV was used. The spin-polarized PBE 56 
generalized-gradient XC functional was employed. An 
un-solvated and hydrogenated gas-phase dimer model 
was extracted 57 from the a (+)Cu(II)Pc 2 polymorph 
structure, with a stacking angle of 65.1 degrees and a dis- 
tance between molecular planes of 3.42 A, giving a lateral 
offset of 1.58 A, as reported from transmission electron 
diffraction analysis described in Ref. 58. A simulation 
cell of 30Ax3oAx20A provided an interatomic spacing 
between periodic images of at least 13.5 A in plane and 
16.5 A out of plane. 



A. Magnetic dipole moments 

The open-shell singlet fragments of the Cu(II)Pc2 sys- 
tem consist of single spins, i.e., a moment of 1 [Lb on 
each copper centre, aligned antiparallel with respect to 
each other. Since, however, approximate XC functionals 
may lower the energy by delocalizing and partially oc- 
cupying magnetization-carrying orbitals 43 , a diminished 



value for the local moment is often recovered in prac- 
tice. The DFT+J7 method seeks to ameliorate this con- 
dition in two complementary ways, that is by introducing 
a derivative-discontinuity to the energy functional which 
penalises fractional occupancies of the spin-orbitals de- 
fined by the subspace projections and also, in doing so, 
by effectively constraining the Kohn-Sham spin-orbitals 
to more closely resemble the (usually more localized) spa- 
tial form of the correlated subspace. 

In spite of this, the correlated subspace projected mag- 
netic dipole moments, shown in Fig. 7, indicate that the 
DFT+C/ method does not effectively localize the magne- 
tization density to the copper 3d manifold for any rea- 
sonable value of the U parameter. Using conventional 
hydrogenic Hubbard projectors, with our best guess for 
the radial profile 40 , we see that there is only a very slight 
increase in the local moment with U. Switching to self- 
consistent projectors in the "tensorial" representation, 
we find that the moment is effectively [/-independent and 
reduced with respect to the hydrogenic result. 

Conversely, the "dual" representation yields a greater 
magnetic moment than the "tensorial" representation, by 
approximately 0.1 /is at [7 = eV, increasing steadily 
at a rate of ps 0.02 fjtBeV . The reason for this discrep- 
ancy, and failure of DFT+U in this regard, is understood 
via the atom-decomposed Mulliken analysis of the mag- 
netization density, which gives 0.10 — 0.12 fig on each ni- 
trogen atom which is a nearest-neighbour to copper, irre- 
spective of either the representation or the U parameter. 
Notwithstanding their adaptation to the molecular en- 
vironment, the self-consistent NGWF projectors remain 
predominantly on the copper ion and do not have suf- 
ficient weight on the neighbouring in-plane nitrogen 2p 
orbitals to capture the magnetization density associated 
with them. As a result, in the same manner as the con- 
ventional projectors, they fail to retrieve the magnetiza- 
tion to the copper 3d x 2_ y 2 orbital within DFT+C7. The 
"dual" representation, however, partially overcomes this 
obstacle, due to the dual Hubbard projectors extending 
over all of the delocalized states in the system including 
those that contribute to the magnetization density. 



B. Kohn-Sham eigenstates 

The accepted understanding 51 ' 54,55 of the spectro- 
scopic nature of the gap in the copper phthalocyaninc 
monomer is that the HOMO level is dominated by a 
doubly-occupied a\ u orbital which consists of a super- 
position of Carbon p z orbitals delocalized on the pyr- 
role rings of both monomer units, while the spectroscop- 
ically correct LUMO level is also a delocalized doubly- 
degenerate orbital, of e g symmetry composed of a su- 
perposition of 7r orbitals on pairs of macrocycle Carbon 
atoms. We may expect some minor differences in the 
spectroscopic properties in the dimer system with respect 
to the monomer, due to cr-bonding between moieties, but 
for the main features to be preserved. 



12 





0.64 r 


PS 
3. 


0.62 - 
0.60 


c 

1) 

s 


0.58 - 


o 


0.56 - 


S 


0.54 - 


1) 


0.52 » 


a 

so 


0.50 


ace m 


0.48 - 
0.46 


Q. 

<j-j 


0.44 


P 

C/5 


0.42 - 
0.40 L 



Hydrogenic 


+ 


Tensorial 




Dual 






> 



2 4 6 

Hubbard U parameter (eV) 



FIG. 7. (Color online) The average magnitude of the projec- 
tion of the magnetic dipole onto the correlated subspaces of 
Cu(II)Pc2, plotted as a function of U for various definitions 
of the subspace projection. 



It has previously been shown that, due to self- 
interaction errors, LDA and GGA-type XC function- 
als do not correctly reproduce the qualitative ordering 
of states close to the Fermi level in the monomer 52,54 . 
The DFT+U insulating gap of the dimer system, within 
various representations, is shown in Fig. 8, along with 
the CZ-dependence of the states nearest the Fermi en- 
ergy. For the spin-polarized PBE functional, we find a 
gap of 0.7 eV for the dimer, whose nature is a charge- 
transfer excitation between b\ g orbitals on either moi- 
ety. The b\g orbital is that which carries the magnetiza- 
tion density in the dimer, consisting primarily of copper 
3d x 2_ y 2 er-bonded to in-plane nitrogen 2p. The repre- 
sentation dependence of the HOMO-LUMO gap follows 
the same trend as the local magnetic moment, due to the 
DFT+t7 correction to the Coulomb-repulsion gap being 
somewhat augmented by an associated enhancement to 
the exchange splitting. In the case of the HOMO or- 
bital, a small value of U is needed to push the singly- 
occupied big state to its spectroscopically correct posi- 
tion below the ai u state and the effect is rather more 
strongly pronounced in the "dual" representation than 
in the spatially-localized methods. 

The incomparability between Kohn-Sham eigenener- 
gics with cither experimental optical or photoemission 
spectra notwithstanding, it is perhaps worth noting some 
similarities and differences between our computed Kohn- 
Sham levels for the dimer system and the gas-phase ul- 
traviolet photoelectron spectra (UPS) and x-ray absorp- 
tion near-edge structure (XANES) reported in Ref. 55. 
Turning first to the valence band edge, the UPS spectra 
confirm that single-molecule CuPc possesses a doubly- 
occupied HOMO of pyrrole-delocalised a\ u character, 
while a further weak ionisation peak at 800 meV above 
that is consistent with the magnetisation-carrying cop- 
per 3d a; 2_ !/ 2-based bi g orbital seen, albeit substantially 
closer to the valence band edge for moderate U values, in 
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FIG. 8. (Color online) The HOMO-LUMO energy gap (top) 
and the energy levels adjacent to the Fermi energy (bottom) 
of Cu(II)Pc2, plotted as a function of U. In the bottom 
panel, solid lines show energy levels of states of predominantly 
Cu-centered bi g character, and to which the DFT+17 correc- 
tion strongly applies, while dashed lines show energy levels of 
states of more delocalized nature. 



Fig. 8. This may suggest a suppression of spin-splitting 
between the bi g levels in the dimer over the monomer sys- 
tem, or may be due to an underestimation of experimen- 
tal transition energies which is not sufficiently alleviated 
by DFT+?7. In the case of the lowest conduction bands, 
carbon K edge XANES spectra assign a large pyrrole car- 
bon character to the LUMO, consistent with a delocalised 
and degenerate e g type orbital. Due to the differing lo- 
calisation regions of the carbon Is and singly-occupied 
bi g orbital, as noted in Ref. 55, such excitations are ne- 
glected and so we cannot compare our prediction that 
the latter orbital lies somewhat below the e g level for all 
except the "dual" representation at high U values. We 
agree on the proximity of these latter levels with previous 
monomer calculations using the PBE functional 53,54 . 

The "tensorial" and "hydrogenic" representations have 
similar effects, as expected; the effect of projector self- 
consistency is rather small in this system. In the case 
of the virtual orbitals, the localized bi g character of 
the LUMO persists for the "tensorial" and "hydrogenic" 
methods, which agree quite closely, while U > 6 eV is 
sufficient to expose a delocalized e g orbital as LUMO 
in the "dual" representation. There is necessarily some 
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small perturbative effect on dclocalized orbitals induced 
by changes to those which are DFT+[/ corrected, which 
is evident in all projection techniques but noticeably 
stronger in the "dual" representation. 

The overall result is that for this hybridized correlated 
system, the "dual" representation recovers the expected 
magnetic dipole moment with significantly more success 
than the fully localized projections. The spectroscopic 
nature of the insulating gap is also recovered to a greater 
degree for a given value of U. We would contend, how- 
ever, that it does so for reasons not expected in the 
DFT+U method. Specifically, where the local magnetic 
moment as measured by the "dual" projectors increases 
with increasing U, the spatial distribution of this increase 
is made up both of the region immediately surround- 
ing the copper ion and spatially diffuse contributions, as 
opposed to the "tensorial" or conventional orthonormal 
"hydrogenic" contributions, with which we are guaran- 
teed to include only subspace-localized densities. 



VI. CONCLUDING REMARKS 

We have presented a revised formalism for the con- 
struction of projection operators, and consequently the 
occupancy matrices, of strongly-correlated subspaces us- 
ing nonorthogonal Hubbard projector functions in ab 
initio methods such as DFT+U and DFT+DMFT. In 
contrast to the previously proposed "full" 13 , "on-site" 14 
and "dual" 8 representations, our "tensorial" definition 
preserves the important property of tensorial invariancc 
in the total occupancy of each subspace, the total en- 
ergy and the ionic forces, by construction. The required 
expressions for the tensorially- invariant DFT+£/ energy 
functional and the resulting potential and ionic forces, 
have been presented. 

Localized nonorthogonal basis functions for Kohn- 
Sham states are frequently used to represent the Hub- 
bard projectors, in practice, either for reasons of com- 
putational convenience or to achieve projector self- 
consistency 15 . We have shown, however, that it is in- 
appropriate to continue to identify the dual space (and 
the metric tensor) of the basis functions with the dual 
space of Hubbard projectors on each site. For molec- 
ular systems, in particular, the unexpected discrepancy 
with respect to orthonormal projectors that is thereby 
introduced may be significant. The resulting projector 
duals (contravariant vectors) are unsuited to construct- 
ing a correction for localized correlation effects, generally 
being delocalized across the entire simulation cell. When 
using delocalized projector duals, moreover, a tensor- 
incompatible symmetrization of the projection operator 
is needed to ensure a Hermitian potential. This may re- 
sult in unphysical occupancy matrix elements and an un- 
controlled action of the corrective potential which it de- 
fines. Put simply, additional non-local corrections are in- 
troduced in the "dual" representation which are extrane- 
ous to the requirement of accounting for the nonorthog- 



onality of the Hubbard projectors. 

Our tensorial formalism may be implemented in any 
methodology which makes use of a nonorthogonal set 
of functions to define each correlated subspace. Since 
it inherently preserves the spatial localization of Hub- 
bard projector duals, it is also less computationally ex- 
pensive and simpler to implement in linear-scaling meth- 
ods, in practice, than the "on-site" or "dual" represen- 
tations which employ delocalized dual projectors. In or- 
der to alleviate the remaining arbitrariness in DFT+U 
and related methods in the nonorthogonal case, the ten- 
sorial formalism may be combined with both a projec- 
tor self-consistency algorithm 15 or any one of a number 
of available first-principles methods for the U parame- 
dians, 21-24. latter remains as an avenue for future 
investigation. 

It is our hope that we have dispelled some of the am- 
biguities surrounding this topic which we feel have arisen 
inevitably as a result of the neglect of the invaluable ten- 
sor notation. As the use of linear-scaling ab initio ap- 
proaches becomes increasingly widespread, we envisage 
that this work may aid the routine implementation of 
sophisticated functionality in the nonorthogonal bases, 
obviating the expenditure of explicit orthonormalization. 
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Appendix A: Orthonormal Hubbard projectors 

Orthonormal sets of Hubbard projectors, as well as 
nonorthogonal sets, may provide a compact and accurate 
representation of the correlated subspaces and we would 
not wish to detract from their value and ease of use. In 
the orthonormal case, the Hubbard projectors equal their 
own duals with respect to their subspace, and the metric 
tensors reduce to Kronecker delta functions. 

If one performs an inverse Lowdin transform 16 from an 
orthonormal set of projectors to a nonorthogonal frame 
using the matrix square root of covariant and contravari- 
ant metrics on a particular correlated subspace, and 
0~ 5 , respectively, then the pre- multiplicative scalar U 
parameter for that site remains identically the same, 
since for each site (if n and n' index orthonormal projec- 
tors and m and m! index their nonorthogonal counter- 
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parts) we have, supposing n, n', m, ml G {1, . . . , M^}, 

^ ^ Unn'nn' — ^ ^ U nn > n " n '" $nn" $n'n'" 
nn' nn'n"n"' 

= E E ^»w» E oLo"^""oi m ,o-^'«'" 

nn' n"n"' mm' 

= E m ' = ^ (/)2 ^- 

mm' 

Thus, when the Coulomb interaction is approximated 
by a pre-multiplicative scalar U times the identity, we 
retain its usual interpretation as the averaged screened 
Coulomb repulsion between densities in the subspace de- 
scribed by the Hubbard projectors, regardless of whether 
or not the Hubbard projectors are orthonormal. 



Appendix B: Invariance under generalized Lowdin 
transforms 

As suggested in Rcf. 9, generalized definitions of the 
Lowdin transform may be envisaged whereby the metric 
tensor is raised to an arbitrary power A, as is its inverse, 
and the canonical Lowdin transform A=\ has the status 
of a special case. Since, however, by construction 

s = V o {I)[A) o {I){ - A)rnn ' = V s (A) s ( - Ahn ' 

n / j nm ^^^^j n^ 

the fully averaged scalar U is invariant under such trans- 
formations, independent of the exponent A, regardless of 
whether the subspace metric tensor O,, or, in the "dual" 
representation case, the metric S 1 ,, on the space spanned 
by all basis functions, S, is used. 

In the latter case of S„ , the generalized Lowdin trans- 
formation exponent A varies the nonorthogonality of the 
representation of the occupancy matrices or, equivalcntly, 
(since the basis set metric S introduces spurious contri- 
butions to the occupancy matrix from across the simula- 
tion cell) the degree of non-locality of the correction. The 



dependence of computed ground-state properties and of 
the Kohn-Sham gap of a variety of materials on A, as 
reported in Ref. 9, demonstrates, in our view, the am- 
biguity of population analysis measures, and hence cor- 
rections such as DFT+J7, which are built from tenso- 
rially inconsistent (necessarily symmetrized) occupancy 
matrices where delocalized Hubbard projector duals of 
the form \ip^ m ) = £ gsl^)^ - ^ 71 ™ are used in the 
construction of the Hubbard projectors. 

The observation of Ref. 9 that the A parameter bears 
influence on computed properties accords well with our 
arguments on the unsuitability of the metric S" (that 
for the basis functions in the entire simulation cell) in 
constructing localized self-interaction corrections such as 
DFT+l/, since that parameter effectively controls the su- 
perfluous spatial derealization of the Hubbard projector 
duals and hence the severity of the tensorial inconsistency 
in the DFT+U functional. By varying A, the occupancy 
matrix for the "dual" representation subject to a gener- 
alized Lowdin transformation, given by 

7,<5eS 

picks up differing non-local contributions (densities from 
outside the correlated subspace). Spurious non-local con- 
tributions are incorporated for all values of A, moreover. 

On the contrary, in the "tensorial" representation, the 
generalized Lowdin transformed occupancy matrix 

oww-"(^|^W)oW( 1 -^'"" i ' ) 

m",m'"eC< 7 ) 

contains no contributions from outside the correlated 
subspace it is describing, for any value of A. Both the 
trace of this matrix and the trace of its square are entirely 
independent of A, since 0^ 1 - A ^ m "' m O^ A ^ mm " = 
(i)m'"m'\ ThuSj by cons truction, the DFT+[7 correc- 
tion is invariant under generalized Lowdin transforma- 
tions and so is unambiguously defined, for a given choice 
of projectors, when the appropriate subspace-local metric 
tensor O" is used to build the projection operator. 
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